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ON EFFICIENT CONSTRUCTION OF STOCHASTIC MOMENT 

MATRICES 

KARRI HAKULA AND MATTI LEINONEN 


Abstract. We consider the construction of the stochastic moment matrices 
that appear in the typical elliptic diffusion problem considered in the setting of 
stochastic Galerkin finite element method (sGFEM). Algorithms for the effi¬ 
cient construction of the stochastic moment matrices are presented for certain 
combinations of afHne/non-afBne diffusion coefficients and multivariate poly¬ 
nomial spaces. We report the performance of various standard polynomial 
spaces for three different non-afHne diffusion coefficients in a one-dimensional 
spatial setting and compare observed Legendre coefficient convergence rates to 
theoretical results. 


1. Introduction 

Continuing progress of efficient deterministic numerical solvers for parametric 
partial differential equation models requires advances in several areas, including, 
the parametrization of the uncertain and spatially inhomogeneous input, the math¬ 
ematical analysis of parametric differential equations, and the actual development 
of the numerical solvers. Advanced engineering applications are emerging as well. 
For instance, in [20] it is noted that more realistic and flexible models for stochas¬ 
tic input (SI) are needed. This work extends computationally efficient modeling 
options in the context of stochastic Galerkin finite element method (sGFEM) with 
a tensor basis, focusing on the specific problem of construction of the associated 
moment matrices. The structure of the induced graph of these matrices, and thus 
computational complexity both in storage and time, depends on the chosen model 
of stochastic input. This problem has been studied in the special case of affine 
stochastic input 0113 E]. In this paper, the construction algorithms are extended 
to cover a much larger class of SI, the non-affine models. These ideas were moti¬ 
vated by the conductivity reconstruction problems where the non-affine stochastic 
input arises naturally via projection of the probability density function of the trun¬ 
cated multinormal distribution onto multivariate Legendre polynomial basis m- 
Under very general assumptions, the computational complexity of the generalized 
algorithms is asymptotically optimal. The numerical results further confirm the 
correctness of the resulting systems. 

Let us introduce the model problem and state the main result more precisely. 
We consider the following model elliptic diffusion problem: Let (fl, E,P) be a 
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probability space. Find random field u G Lp{Q, Hq{D)) such that 


r — V • (a(a;,x)Vu(w,x)) =/(x) in D, 

[ m(w,x) = 0 on dD 


holds P-almost surely for a load / G L‘^{D) and a given strictly positive diffusion 
coefficient a G L°°(Q, x D) with lower and upper bounds aminjamax such that 


(1.2) P ( w G n : 0 < amin < ess inf a(a;, x) < esssupa(a;,x) < Umax I = 1 - 

V xGl> / 

Moreover, the physical domain P is a bounded domain in where d G N, with a 
smooth enough boundary and Pp(n, Hq{D)) is a Bochner space; see, e.g., [18] and 
the references therein for more information. The stochastic input, i.e., the diffusion 
coefficient, in the model problem is usually assumed to be given in the following 
affine form: 

Definition 1.1 (Affine diffusion coefficient). The affine diffusion coefficient is de¬ 
fined as 

(1.3) a(a;,x) = ao(x) + ^ am(x)y„(w), 

m>l 

where {am}jTi>o are some suitable spatial functions and {Ym}m>i is a family of 
random variables. This affine form is often obtained by using the Karhunen-Loeve 
expansion; see, e.g., |2]. 

As noted above, in practical applications we might encounter non-affine diffusion 
coefficients, and for example, in |13j it is stated that many practically relevant 
parametric PDFs are nonlinear and depend on the parameters {l^}m>i in a non- 
affine manner. Formally, we define the relevant concepts as follows: 

Definition 1.2 (Support of a multi-index). The support of a multi-index rj G Nff 
is defined as supp(77) = {m G N | 0}. 

Definition 1.3 (Finitely supported multi-indices). The set of finitely supported 
multi-indices (N)(°)c is defined by 

(Ng°)c = {77 G I |supp(77)| < 00} C 

where we use |supp(? 7 )| to denote the cardinality of supp)??). We define the set 
(Z°°)c in a similar way. 

Definition 1.4 (Non-affine diffusion coefficient). The non-affine diffusion coeffi¬ 
cient is given by 

(1.4) a(w,x) = ao(x)-P ^ a^(x)y^(w), 

mgh 

where S is a subset of finitely supported multi-indices, oq and are some 

suitable spatial functions, and 

00 

y^(a.) := n = n 

m—1 mGsupp(/i) 
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For example, the (typical) instance 

a(uj,x) = ao(x) + ^ am(x)Ym(uj) 


I m>l 


considered in m fits the non-afhne form ( |1.4[ ). In this paper, we refer with non- 
affine to the form (1.4) and we do not consider other non-afhne type diffusion 
coefficients, e.g., log-normal, unless explicitly stated. Notice that the non-affine 
case includes the affine case; by selecting E = {fi G (Ng°)c | /r = e„ for some n G N}, 
where e„ denotes the nth Euclidean basis vector of , the non-affine case reduces 
to the affine case. Here, as in the following, we have used the convention 0° = 1. 


A standard approach is to transform the model problem (1.1) into a parametric 
deterministic form by hrst assuming the family H —>■ M to be mutually 

independent with ranges F^ and associating each with a complete probability 
space {flm, Ym, Pm), where the probability measure Pm admits a probability density 
function pm ■ F^ —> [0,oo) such that 

dEm(cc) — Pmiym)^ym, ym ^ Ym, 

and the a-algebra Ym is assumed to be a subset of the Borel sets of F^. Finally, it 
is assumed that the stochastic input data is finite, i.e., in the affine case we assume 
that there exists M < oo such that Ym = 0, when m > M, which essentially 


truncates the series in (1.3) after M terms, and in the non-affine case (1.4) we 


assume that the set S has a finite number of elements and M as in the affine case 
exists. 

After these transformations, the parametric deterministic weak formulation of 


(1.1) is to find u G Lp(r, Hq{D)) that satisfies 


(1.5) / / a(y,x)Vu(y,x) • Vu(y,x)p(y)dxdy = / / /(x)u(y, x)p(y)dxdy 

JrJ d JrJ D 

for all V G where y := {yi,y 2 , ■ ■ ■) GT := Fi XF 2 X ... and p(y)dy := 

Y\m>i Pm)yrn)^ym- The Unique solvability of ( |1.5| ) follows from the Lax-Milgram 
lemma. 


By discretization, the weak formulation (1.5) reduces to the following linear 


system of equations in the affine case (1.3): 

M 

( 1 - 6 ) I Go <Yi Aq + Gm <S> Am I U = go<Y> fo, 


where {Am}m=o standard finite element (FEM) matrices, /o is a standard load 
vector, {Gm}m>o are stochastic moment matrices, and go is a stochastic load vector; 
see m for details. 

For the stochastic load vector and moment matrices, we recall that each multi¬ 
index 77 G (Ng°)c determines a multivariate polynomial 4?^(y). 

Definition 1.5 (Multivariate polynomial). Let 77 G (Nff’)c. The multivariate poly¬ 
nomial $^(y), also called “chaos polynomial”, is defined as 

00 

®r,(y) = n ^Urr,{ym) = n ^Urn.iym), 

m—1 mGsupp(77) 


(1.7) 
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where {<pm}m>o is a suitably orthonormalized (univariate) polynomial sequence 
with (f>Q = 1 . 

In this paper, we assume that {</>m}m>o is either the Legendre or the (prob- 
abilists’) Hermite polynomial family depending on the family {Lm}m>i as those 
are typical choices in the sGFEM setting; see Appendix [X| for the definition of the 
Legendre and Hermite polynomials. We call the different cases as the Legendre 
case and the Hermite case, respectively. 


Remark 1.6. We consider the Hermite case mainly for completeness as the condi¬ 


tion (1.2) combined with the non-afhne form (1.4) has to be relaxed for the use 


of normally distributed random variables {Ym}m>i', we refer to m 0 na [ig for 
possible relaxations of condition (1.2). We content ourselves by providing the for¬ 


mulas for the stochastic moment matrices in the Hermite case, and in the numerical 
experiments of Section 0 we consider only the Legendre case. 

In the affine case, the stochastic load vector is given by 


( 1 . 8 ) 


■■= J^^aiy)p{y)dy 


and the elements of the stochastic moment matrices are 


(1.9) 


[^o] a ./3 ■= ^a(y)®/3(y)p(y)dy and 

•= J ym‘^a{y)‘^ 0 {y)p{y)dy, m > l. 


where a,l3 G A C (Njf’)c. The index set A is finite and consists of suitably cho¬ 
sen finitely supported multi-indices; see |18j and the references therein for more 
information. In an ideal case, the index set A is chosen so that 


msa 


u/.(x)«>^(y), 


where :=L u(y, •)$^(y)p(y)dy S Hq{D), gives a good approximation to the 
solution of ( |1.5| ); in practice, the index set is often chosen more or less arbitrarily. 
It is known that the stochastic moment matrices {Gm}m>i exhibit a nontrivial 
sparsity pattern |llj . (Gq is a diagonal matrix.) 


Efficient construction of the stochastic moment matrices (1.9) is crucial in the 


sGEEM setting as it easily becomes one of the bottlenecks for the method when 
the number of multi-indices in A increases. One solution is to build the index set 
A so that an efficient construction of the moment matrices is possible; Bieri et al. 
have given such an algorithm in 0 for certain index set types. 


In the non-affine case (1.4), the resulting linear system is (cf. (1.6)) 


( 1 . 10 ) 


Go 0 Aq -k G^ 0 u = go (g) /o. 




where Aq and {A^j^gs are standard FEM matrices and the stochastic moment 
matrices are given by 


( 1 . 11 ) 


[G^ 


n I ®a(y)®/3(y)p(y)dy 


V, 771=1 
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with a and /3 as above. With this notation, the matrix Gm, where m € N, is given 
by G®"* and Gq is given with /i = (0,0,...). As above, these moment matrices 
employ a nontrivial sparsity pattern. 

The efficient construction of the matrices {G'^j^gs is inspired by [HI Algo¬ 
rithm 4.13] for {Gm}m>o- We present a non-recursive version of the original and 
generalize it for {G'^j^gs. This leads to a novel setup for updating the neces¬ 
sary auxiliary information in contrast to [H] (and its research report version [TO])- 
This allows us to form the stochastic moment matrices for the model problem (1.1) 
efficiently when the stochastic input is given in a non-affine form (1.4). 

The computational complexity depends on the set of multi-indices A and the cho¬ 
sen model of the stochastic input. However, the latter part has finite complexity in 
the sense that it has to be fully resolved and thus it is not asymptotic. Therefore, 
the computational complexity is ultimately determined by the set of multi-indices 
only. However, the situation is similar to all sparse versus full considerations. For 
a sufficiently small set of multi-indices it is possible that the stochastic input de¬ 
termines the observed complexity. Asymptotically, we get the complexity estimate 
0(|A|^+®), where 0 < e < 1, and in many practical cases e <C 1. 

The rest of the paper is organized as follows. In Section we give definitions 
for the multi-index sets A considered in this paper and fix some terminology and 
conventions related to multi-indices. Section [^ gives steps for constructing the 
stochastic moment matrices {G^j^gs efficiently. Algorithms for the construction 
of the multi-indices and neighbour matrices required for {G^j^gs are presented 
in Section [^ and our numerical experiments are presented in Section In the 
numerical experiments, we compare the performance of different multi-index sets 
and compare observed Legendre coefficient convergence rates to theoretical results. 
Finally, we conclude with discussion in Section]^ In Appendix [A] we recall some 
essential properties of the Legendre and Hermite polynomials which were employed 
in Sectionj^ and Appendices jBHFjcontain the actual pseudocodes for the algorithms 
presented in Section 


2. Preliminaries 

To improve the readability of this paper, we start by introducing the six multi¬ 
index sets considered in this paper: 

1. isotropic tensor product (isoTP) index set, 

2. isotropic total degree (isoTD) index set, 

3. anisotropic tensor product (aTP) index set, 

4. anisotropic total degree (aTD) index set, 

5. threshold sequence (TS) index set, and 

6. the best M-terms (bestM) index set. 

Standard choices for the index set A are the isoTP and isoTD index sets. Se¬ 
lecting the optimal index set a priori in the sGFEM setting is an active research 
question; see, e.g., [71IH1I5] and the references therein for more information. 

Definition 2.1 (Isotropic tensor product index set). Let N gN and K G Ng- The 
isoTP index set is given by 

a G (N)]°)c I max < AT; = 0, n > TV 



( 2 . 1 ) 


isoTP(A^, AT) 
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Definition 2.2 (Isotropic total degree index set). Let N G N and K G Ng- The 
isoTD index set is defined as 


( 2 . 2 ) 


isoTD(7V, itT) = {a G (Ng°)c 


N 

E 

n—l 


an < K;an = 0,n > N 


The isotropic index sets are often preferred due to their easiness of construction. 
However, the cardinalities of the isoTP and isoTD index sets are {K + 1)^ and 
i^K^) ’ respectively, and hence the isotropic index sets grow rapidly when either N 
or K is increased. Therefore, as the isoTD index set grows slower than the isoTP 
index set, the isoTD index set is often preferred over the isoTP index set. 

The isotropic index sets give equal weight to the N dimensions, and as in some 
applications we might want to prefer some of the dimensions, the anisotropic ver¬ 
sions of the isotropic index sets are also used. 

Definition 2.3 (Anisotropic tensor product index set). Let N G fi, K G Ng, and 


p. = {pi,...,pn) G 


be a vector of positive weights. The anisotropic tensor 


product (aTP) index set is given by 


(2.3) aTP(iV, AT, ^) = <^ a e (N“)c I max /x„a„ </TminAT; «„ = 0, n > iV 

Here and in the following, IR+ := {a; G M | a; > 0} and Pmin ■= min(/i). Similarly, 
we define ya-max = max(/i). 

Definition 2.4 (Anisotropic total degree index set). Let N G N, K G Ng, and 
p = {pi ,..., pn) G be a vector of positive weights. The anisotropic total degree 
(aTD) index set is defined as 

(2.4) aTD(A^, A:,^) = |a e (N“)c I ^ < Amin AT; a„ = 0, n > ivj . 

The weight vector p is used to give different weights to the N dimensions, and 
thus the cardinalities of the anisotropic index sets depend heavily on the chosen 
weight vector. Without loss of generality, we assume that the weight vector p used 
in the aTP and aTD index sets is non-decreasing, i.e., Amin = Ai ^ A2 < • • • < Aw- 
(If it is not the case, we permute the N dimensions according to the weight vector.) 
We refer with TP and TD to both anisotropic and isotropic versions of the tensor 
product and total degree index sets, respectively. Notice, that by using a constant 
weight vector, the anisotropic index sets reduce to the corresponding isotropic index 
sets. 

Closely related to the TD spaces is, what we call, the threshold sequence (TS) 
index set. 

Definition 2.5 (Threshold sequence index set). Let 0 < e: < 1 and a be a non¬ 
increasing sequence of nonnegative real numbers bounded above by one and tending 
to zero, i.e., 

A G cg := < A = (A1) A2, • ■ •) G | lim p^ = Q and l> pi> P2> ■ ■ 

L m—¥oo 

where Kg := {a; G M | a; > 0}. The TS index set is given by 


}• 


(2.5) 


TS(A,e) = laG{N^)c 


a“ := n 

m^N 


n G 

'm.Gsupp(a) 


> £ 
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In [9], essentially the same index set as TS(^,c) is denoted with Ae(/i). (In [9], 
is replaced with K^.) 


By noticing that the aTD set condition 

N 


^ ^ ^ Mmin-A" 

n—1 


can be written as 


by defining 


M = 



N 




> e 



and e = exp{—K), 


an aTD(IV, K, /i) index set can be considered as a TS(/i, e) index set by permuting and 
padding p, so that the vector is non-increasing and belongs to Kg®. In most cases, 
we are interested to construct an index set with less than some maximum number 
of multi-indices and not interested in the actual parameters used to construct the 
index set. Hence, in the anisotropic cases it is often enough to define the weight 
vector /i up to a positive constant ^ 7 ^, where c G K+, by using the vector 


9 = 


. Ml 

/^mir 


Mmin 


i^. 


For example, the aTD index set condition (2.4) can be written as 

exp - ^ gnan > exp{-cK) 


n—1 


and by considering exp(—ciF) as a real number between zero and one independent 
of c and K, we can control the number of elements in the index set by varying the 
value of exp(—ciF). To summarize: 

• The aTD index set can be considered as the TS index set. 

• The TS index set allows us in many practical cases to define the weight 
vector p in the TD index set up to a constant. 

Similarly in the aTP case, it is often enough to define the weight vector /x up to a 
positive constant. 

Finally, the bestM index set is obtained from an existing index set. 


Definition 2.6 (Best M-terms index set). The bestM index set is obtained from an 
existing index set A by associating each element of A with a real valued coefficient, 
usually the norm of after which the elements of A are sorted in a non-increasing 
order based on the corresponding coefficients. The first M-indices from the resulting 
ordered sequence form the bestM index set. 


Using the bestM index set in practice is often impractical as to construct the 
index set, i.e., to obtain the norms of we first compute an sGFEM solution 
corresponding to a large index set and there is often no benefit in considering 
smaller index sets than the large index set already used to construct the sGFEM 
solution. Hence, we focus on the other index sets when discussing the obtained 
results and use the bestM index set only as a reference index set. However, there 
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are at least two cases (not considered in this paper) where it might be useful to 
actually consider the bestM index sets. First, evaluating the sGFEM solution might 
be expensive and one option to reduce the evaluation cost is to reduce the size of 
the index set using the bestM index set. This technique will, in general, increase 
the approximation error but by dropping the multi-indices corresponding to the 
smallest u^, we expect the approximation error to increase only slightly. Second, 
if we are using an adaptive scheme, like reservoir sampling [TU], to construct the 
index set, it may be useful to consider bestM index sets; constructing the index set 
using randomized algorithms is left for future studies. 

Let us close this section with some terminology and conventions related to multi¬ 
indices. In this paper, the addition, subtraction, and product of multi-indices is 
carried out componentwise and, if necessary, the elements of (Nf(°)c are considered 
as elements of (2°°)^. Moreover, we use the following definitions for the length and 
parenthood condition of a multi-index. 

Definition 2.7 (Length of a multi-index). Let a G (N^)c- The length of a is 
defined as 


length(a) = max(supp(a) U {0}), 

where the support of a multi-index was defined in Definition |1. 2 1 

Definition 2.8 (Parenthood condition of a multi-index). Let a and /3 be two multi¬ 
indices. We call a the parent of /3 or /3 the child of a if there exists n S N greater 
or equal to length(Q!) such that /3m = Om + i^mn for all m G N, where Smn is the 
Kronecker’s delta. 

In the algorithms presented in this paper, we store the parenthood relations 
between the multi-indices of a index set A in a variable denoted with PtC in such 
a way that the children of a S A are given in the vector PtC[a]. (Actually, the 
children of a are given in the vector PtC[a], where a G N is an index corresponding 
to a, but for notational reasons we write PtC[Q;] in the text and use the correct 
form PtC[a] in the actual pseudocodes.) Moreover, each multi-index a G (Ng°)c or 
vector a G (Z°°)c is stored in the presented pseudocodes in a sparse format 

a = {(m,am) | ^m 7^ 0} 
and each multi-index a G TS(/i,e) is stored as 

= {{ (m, Otra) I O^m 7^ 0}; ? 

for example, the multi-index (0,1, 2,0,...) would be stored in TS(/x, e) as 

{{(2,l),(3,2)},^i/r2} 

if it belongs to the set, i.e., if ^ In the pseudocodes, we denote with 
a[TO], where m G N, the mth pair in the sparse format of a. Hence, in general, 
a[m] does not necessarily correspond to am- Moreover, we denote with a[—1] the 
last pair in the sparse format of a, and we use short-circuit evaluation of Boolean 
operators, i.e., the second argument of a Boolean expression is evaluated only if the 
first argument does not suffice to determine the value of the expression. 
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3 . Construction of stochastic moment matrices 


Combining Definition 1.5 to (1.111, we can write the element of as 


[Gf‘ 


I Q ;,/3 


n 2^™"* ] ®a(y)®/3(y)p(y)dy 


MS’'”] 

n 

mGsupp(/j,+Q:+/ 3 ) 


(n j f n 

\m—l / \m—l 


Aym)\ p(y)dy 


where a and /3 are as before and the matrix is given by (here, am + 1 and 
fim + 1 are used as general indices and are independent of m) 


[^m”] + ~ J (y7Ti)<('/3„ (2/m)Pm(l/m)d?/rn- 

We do the following assumptions in the Legendre and Hermite cases, respectively: 
Assumption 3.1. In the Legendre case, we assume that 
Lm = [-1,1] := To and 

for all m > 1. 

Assumption 3.2. In the Hermite case, we assume that 

r„=K:=ro and Pm{y) = exp{-y'^/2)/■■= Po{y) 
for all m> 1. (Recall remark 1.6 close to the end of Section^) 


Pmiy) = 1'-= Po{y) 


Due to Assumptions 3.1 and 3.2 G^ = G„ for all A: e No and to, n > 1 (each of 
the stochastic dimensions has the same range and probability density), and hence 
it is enough to consider only the matrix given by (here again I + 1 and to + 1 
are general indices) 


[k’^ 


■= / y^((iiy)4>m{y)poiy)dy, 

'To 


J Z+l,m+l 

where k,l, and to are nonnegative integers, as we can write 

P'l) [G'‘l«,s= n 

mGsupp(/^+Q:+/3) 

By recalling that we are using orthonormal multivariate polynomials, is an 
identity matrix, and therefore G^ can be computed as 


(3.2) 


n [K>^ 


[G^ I fl = •( mesupp(/j) 


Ja„ + l,/3„ + l 


, when 


J a ,/3 


0 , 


otherwise. 


To use (3.2), we need to locate the nonzero elements in G^, which amounts to 
understanding the structure of as is evident from (3.1). 

Using the triple integral and inversion formulas for the Legendre and Hermite 
polynomials from Appendix we can locate and compute the nonzero entries in 
the matrix K^. 
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Theorem 3.3. In the Legendre case, we have 

k 

(3.3) 
where 




n=0 


In = l{k- 


-n IS eveni 


and 


Proof. 


fk\ , /-- {k — n — l)\\ 

In! \/2n + 1-^- 

\nj (/c + n + 1)!! 

L{l,m,n) = \l{21 + l)(2m + l)(2n + 1) ^ 


0 < n < k, 


I m n 
0 0 0 


[K^ 


I Z+l,m+l 


y^Li{y)L^{y)po{y)dy 


J -1 

k 

n=0 


Ln{y)Li{y)Lni{y)po{y)<^y = 


k 

n=0 


where we have used first Theorem IA.5I and then Theorem IA.41 


□ 


Combining the selection rules of and L{l,m,n) from (A.5| and (A.3), respec¬ 
tively, we know that the summand in (3.3) is nonzero when 

• k — n is even, 

• I + m + n is even, and 

• 1^ — m| < n < Z -|- TO, where 0 < n < k. 

From the first condition, we see that if k is even, n is also even and then the second 
condition imposes Z -|- to to be even as well. Hence, Z and to have the same parity 
when k is even. Similarly, when k is odd, n is also odd, and as a result Z and to 
have the opposite parity. From the last condition, we see that is a symmetric 
{2k -b l)-diagonal matrix. To sum up, is a symmetric {2k + l)-diagonal matrix 
with odd diagonals zero when k is even and even diagonals zero when k is odd. For 
example, the common case W used in the affine case is given by m (correcting a 
misprint) 

' TO = Z -b 1, 


V(2i-1)(2/-Hl)’ 

[W] ={ — 


(2i-3)(2Z-l) ’ 


0 , 


TO = Z — 1, 
otherwise, 


where Z,to G N. Similar results hold for the Hermite case. 
Theorem 3.4. In the Hermite case, we have 

k 


(3.4) 

where 




n—0 


K = hk- 


■n IS even 


^-s/ul (Zc — n — 1)!!, 0<n<k, 


H{l,m,Tl^ I{2g is even} ^{\l—m\<n<l+m} 


I 


m 


9 - nj \g - nj \g - { 


and 
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with g = (I + m + n)/2. 

Proof. Similar to the Legendre case; use Theorems |A.9| and |A.8| 


□ 


By noticing that the selection rules in the summand of Theorem |3.4| are the same 
as in the Legendre case, has the same structure as in the Legendre case. The 
common case is in the Hermite case 


VI, 


m = I + 1, 


^ m = l-l, 


0 , 


otherwise, 


where l,m S N. 

Using the structure of and (3.11, we see that the element ^ i® nonzero 

if the following condition holds for all m > 1: lom —/3m| < Mm and am and Pm have 
the same (opposite) parity when g,m is even (odd). This suggests the following for 
the Legendre and Hermite cases: Associate each multi-index g G (N((°)c with the 
following set of weights 


1U(m) := {aP I a G P{g,) and p G F{V)} C (Z°°)c, 


where 


and 


pv) ■■= n 


2>1 


{!}, Mm = 0 or = 0 for all n < m, 
{1,-1}, otherwise. 


^’(m) := n 

ra>l 


{0,2,-■ ■ ,Hm}, Mm is even, 
{1,3, * * * , Mm }, t^m is odd. 


Now, the condition for the nonzero elements in can be given as: 

(3-5) [G^] „,/3 ^ ^ a,/3 ^ ® a,/3 ^ 

where TV’" is the neighbour matrix given by 


and is the summed matrix 


1, a ± u> = /3, 
0, otherwise. 


:= N^. 

wSiW{fi) 

Finally, the stochastic moment matrices for finitely supported multi-index 

sets S and A can be computed in the following way: 

(1) Construct the neighbour matrices {A^’"}, where w G We := (J W{fi). 

(2) Compute the summed matrices using {A^’"}u,gvus- 

(3) Construct the matrices {AT^j, where fc G < 1,..., max(max(p,)) >, using 


Theorem 
(4) Compute 


3.3 


or Theorem 


3.4 


Ate; 


tEe matrices {G^j^gs using (3.2) and (3.5). 












12 


HARRI HAKULA AND MATTI LEINONEN 


Here, iV’", S^, and are symmetric (sparse) matrices of dimension |A|, and 
is a symmetric (2k + l)-diagonal matrix of dimension max(max(a)) + 1. 

cxGA 

The second, third, and fourth steps from above are straightforward to carry out. 
However, the first step, i.e., constructing the neighbour matrices is a 

nontrivial task. Straightforward methods to construct the neighbour matrices such 
as using the definition (3.6) directly or writing w = ±(a — /3) and looping over a and 
fi result in algorithms where we essentially loop through We once for each {a,/?} 
pair, i.e., the construction requires 0 (|Wh||AP) operations. In the next section, we 
give an algorithm based on m which can be used to construct hierarchical multi¬ 
index sets so that constructing neighbour matrices for these index sets require 
C>(|Wh||A|^+'^) operations, where 0 < e < 1; the exact value of e depends on S 
and A and is much less than one in typical instances. The speedup is obtained by 
constructing the multi-indices in A in a such way that locating a multi-index in A 
does not necessarily require us to loop over entire A. 

The cardinalities of the sets H(^) and F(/i) are given by 

2I"UPP(M)|-1, length(/.) 


IAm)I = 


T 


/r = 0, 


and |F(/r)| = 


m—1 




respectively, from which we get the following estimate for the cardinality of W{fi): 

f 2lsupp(A»)|-i nS'""’ + ij 

1 , 


(3.7) \W(y)\ < \P(MF(k)\ = 


^ ^ 0 , 
^ = 0 . 


From (3.7), we obtain a (conservative) upper bound for the number of (sparse) 
matrix additions required for the computation of the summed matrices in the second 
step: 




\W(^i)\ < 1 


length(/^) 


AiGH 

m/O 


Mm 

lY 


1 


<1 + 2.2 

/iGH 


M-1 ( max(/x) 

V 2 


-f 1 


M 1^1 


< 


H(: 


max(max(/i)) -|- 2 

V 


M 


Notice that the above complexity estimate depends entirely on the stochastic input, 
and hence the time spent on the second step is essentially imposed by the prob¬ 
lem formulation; the simpler the stochastic input, the less time we are spending 
on computing the summed matrices. In the simplest case, i.e., in the affine one, 
IT(m) = {m}i and hence and the second step can be skipped. 

By combining the complexity estimates for the first and second steps, the overall 
complexity of locating nonzero elements in the stochastic moment matrices is 


O |IFh||A| 


l-|-£ 


( max(max(/r)) 

V 


M^ 


If we assume the stochastic input in the problem formulation, i.e., S and M, to be 
fixed, the complexity estimate amounts to O (|A|^+'^). 

Employing the neighbour matrices is profitable if the time used to construct 
them and the summed matrices is less than the time used to evaluate the zero 
elements in the stochastic moment matrices. If the stochastic moment matrices are 
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sparse, the time spent on evaluating the zero elements is essentially (!l(|Ap) using 
(3.1) directly and 0(|A|^+'^), where 0 < e < 1, if we employ the neighbour matrices. 
However, if the neighbour matrices are dense, the time spent on evaluating the 
zero elements using (3.1) directly is practically zero, whereas the time to construct 
the neighbour matrices has not changed, and hence it is better to construct the 
stochastic moment matrices using (3.1) directly. Be that as it may, by employing 
the neighbour matrices, we have been able to reduce the stochastic moment matrix 
construction time from several hours to a few minutes in many practical cases. 


4. Constructing hierarchical multi-indices and neighbour matrices 

In this section, we present algorithms that can be used to construct the TS,TD, 
and TP index sets and the corresponding neighbour matrices {A^“'}u,gvUH • We give 
both conceptual (Algorithms and and pseudocode (Algorithms and im¬ 
plementations of the algorithms. We consider hrst the TS index set and then give 
the small modifications required for the TD and TP index sets. We conclude by 
considering other similar index sets. 

4.1. Threshold sequence index set. The algorithm given in Appendix [Bj i.e., 

Algorithmic which is the stack-based version of the recursive algorithm presented 
in pro], can be used to construct the index set TS(/i, e). Subsequently, the neighbour 
matrices can be constructed using Algorithm |C from Appendix [P] 

Algorithm [C takes a sequence fj, G cq and a tolerance £ > 0 as input and returns 
the index set TS(/r,e) together with the corresponding parenthood information (cf. 
Definition |2.8[ ); the parenthood information is required in the construction of the 
neighbour matrices in Algorithm [C 

4.1.1. High-level description of the TS index set construction. To ease the under¬ 
standing of Algorithm [C a high-level description of the algorithm is given in Algo¬ 
rithm and explained more in detail below. The line(s) indicated at the end of the 
lines in Algorithm correspond to the relevant pseudocode lines in Algorithm [C 

Let /3 be a child of a multi-index a. The multi-index /3 is a feasible child of 
a if the index set condition > e holds. After a has been added to the index 
set under construction, the feasible children such as /3 are added to the index set 
in the decreasing order of length before any other multi-indices are added to the 
index set. Moreover, the parenthood information of the just added multi-index is 
recorded in a natural tree-format. The index set TS(/r, e) is always initialized with 
zero multi-index. 

The feasible children are processed in the decreasing order of length to ensure 
compatibility with the recursive version of the algorithm presented in nni and 
to facilitate efficient neighbour matrix construction later in this section. (Any 
consistent choice here would work with some minor modifications to the algorithms.) 

The algorithm is guaranteed to stop at some point as there is a finite number of 
multi-indices for which the condition > e holds; see [5] for the analysis of the 
running time and for bounds of the size of the index set TS(/r,e). The algorithm 
generates the set TS(^,e) essentially due to the fact that if < e then, by 

recalling that p, S cq, it holds that pL°‘p.m < £ for all m > n as well. Hence, if 
a child of the current multi-index is not feasible, we know that any sequence of 
one-increments to the child would not produce a multi-index that would belong to 
the index set, and therefore there is no need to consider it. An example run of 
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Algorithm with intermediate steps in the spirit of Algorithm [T] is presented in 
Appendix [Cl 


Algorithm 1: High-level description of the TS index set construction. 

1 set current multi-index to (0, 0,...) and add it to TS(/i, e); // 1-5 

2 while True do 

push all feasible children of the current multi-index into a stack in the 
increasing order of length; // 7-12 

if stack is empty then return TS(/r, e) and the corresponding parenthood 
information; // 13 

set the current multi-index and the associated variables to the top state in 
the stack and pop the stack afterwards; // 14-18 

add the current multi-index to TS(/i,e) and record the corresponding 
parenthood information; // 19-23 


7 end 


4.1.2. High-level description of the neighbour matrix construction. Algorithm [^ 
given in Appendix [ d[ can be used to construct the neighbour matrices {A’"}u,gWh- 
The algorithm takes the weights HA, the index set TS(/r,e), and the corresponding 
vector p, real number e, and parenthood relations as input, i.e., the input and 
output of Algorithm]^ and returns the neighbour matrices 

The high-level description of Algorithm is given in Algorithm in a similar 
format as Algorithm and is essentially: for each w G HA loop over each rj G 
TS(^,e) and if 7 = 77 — w belongs to TS(/r,e), locate 7 G 1S{p,e) and add the 
corresponding contribution to A™. This method requires an efficient way for both 
checking whether 7 belongs to TS(/r,£:) and locating 7 in TS{p,e). 

A simple and efficient test for 7 = 77 — w ^ TS(/i, e) is given by: 

7 ^ 13(77, e) 7 m < 0 for some m or p? = 77’'““ = p^/p^ < e. 

We need to check 7 for possible negative values as after subtracting w from rj we 
may end up with negative values in 7 . Moreover, notice that we are storing the 
value p'^ already in the sparse format of 77 in 18(77, e), and hence it is enough to 
compute the value 77“ once for each w G HA. 

The efficient location of 7 in 13(77, e) is based on the tree structure provided 
by the parenthood relations. We can locate a given multi-index 7 G 18(77, e) by 
starting from the root of the tree, i.e., ( 0 , 0 ,...), after which we move down the 
tree by processing the position-value pairs (m, 7 ^) in the sparse format of 7 in the 
increasing position order using the procedure: 

( 1 ) Move one step down the tree to the child with the length m. 

(2) Move 7 m — 1 times down the tree selecting the child with the length m. 

In our implementation, the second step corresponds to selecting the first child in 
the child list of the current multi-index m times. 
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For example, the multi-index (2, 0,1,0,...) = {(1,2), (3,1)} in the example run 
of Algorithm 1^ presented in Appendix is located through the steps: 

_(F2)_ 

( 0 , 0 , 0 )^( 1 , 0 , 0 ) ^( 2 , 0 , 0 ) ^( 2 , 0 , 1 ). 

(3.1) 


Here, we have used the consistent recording order of the parenthood relations and 
the fact that only nonzero indices are stored in the sparse format. 


Algorithm 2: High-level description of the neighbour matrix construction. 

1 for w G Ws do 

// 1-3 

2 

for rj G TS(/r, e) do 

// 4-6 

3 


if r] — w ^ TS(^, e) then continue; 

// 7-8 

4 


locate r] — w in TS(p,,e); 

// 9-22 

5 


add contribution to A*"; 

// 23 

6 

end 


7 end 

8 return {N^}w^Ws\ 

// 26 


4.2. Total degree and tensor product index sets. Algorithms and can 
be readily modified for the TD and TP index sets; in Algorithm we change the 
definition of “feasible children” and in Algorithm the set condition for rj — w. 

In principle, the TD case does not require any modifications as it is possible to 
reduce the TD case to the TS case as explained in SectionHowever, in the isoTD 
case we might want to modify the algorithms slightly to avoid any possible problems 
arising from the floating point arithmetic; we perform a global replace of TS(/r, e) to 
isoTD(iV, AT) and the modifications presented in Appendix to Algorithms an d 
l^to bring the algorithms compatible with the set condition (2.2) instead of (2.5). 

In the TP case, we need to consider only the aTP index set as the aTP index set 
reduces to the isoTP index set when ^ is a constant weight vector. As in the isoTD 
case, we perform a global replace of TS(/r, s) to aTP(iV, AT, /r) and the modifications 
presented in Appendix [F| 


4.3. Complexity estimates. From the high-level descriptions of Algorithmsand 
we see that the workload of the algorithms is C(|A|) and 0 (|IFh||A| max^gA |Q!|), 
respectively, where |a| = "^he value of e in the estimate C>(|Wh||A|^+'^) 

presented at the end of Section can be approximated by equating 

(4.1) max lal = lAI*^ e = log ( max |a| ) / log(|A|). 

aeA yaeA j 

From the estimates 


(4.2) 


max |a| < max |a| = A^AT, 

a6aTP(A,A',/j) a61soTP(A',if) 

max |a| < max |a| = AT, 

aGaTD( Q:GisoTD(A^,iC) 











16 


HARRI HAKULA AND MATTI LEINONEN 


and 


(4.3) 


Mmin 

Mmax 

N + 


K 


N 


+ 1 < I aTP(7V,K,/i)I < I isoTP(7V, K)\ = {K + lY, 


Mmax 

'n 


K 


< I aTD(iV, K,n)\<\ isoTD(iV, K)\ = 


N + K 
N 


we get by combining (4.2) and (4.3) in (4.1) and letting K ^ oo the upper bound 


lim e < — 
K—^oo N 


for the TP and TD index sets. In most cases N = M, and hence the complexity of 
the neighbour matrix construction is approximately O j for the TP and 

TD index sets. 

Timing results for a diffusion coefficient of the form 


(4.4) 


j(y,x) = 1 + 


M 

1 + a„ 

m—1 




and an isoTD(M, AT) index set with different values of M,p, and K are presented 
in Figure The tables on the left of Figure portray, from top to bottom, the 
estimated value of e in the complexity estimate C>(|IFh| lAj^^*^) together with the 
theoretical value M~^, the number of (sparse) matrix additions required for the 
computation of the summed matrices, i.e., I^(m)I) maximum value 

of K considered when estimating e for different M and p. The image on the right of 
Figure[^shows the time required to compute the neighbour matrices for M = 6 and 
p = 2,4, 6 together with the line used to estimate e. The plot markers correspond to 
different values of K, and the timings were carried out with Intel® Xeon® Processor 
E3-1230 V2 (8M Cache, 3.30GHz) with 15.6 GiB of memory. 

The obtained values of e depicted in Figure are in agreement with the theoret¬ 
ical limits. Increasing p increases the y-intercepts of the fitted lines in the timing 
plot which is in agreement with the complexity estimate as increasing p amounts 
to a larger HG- 


4.4. Other index set types. It is possible to consider other index sets, say A, than 
the TS, TD, and TP index sets using Algorithms]^ andwith (minor) modifications. 
What is required is a method to check whether a given multi-index is an element 
of A and some structure like the parenthood relations to locate it. Gurrently, the 
tree structure provided by the parenthood relations requires that the index set A 
can be constructed by adding an initial multi-index to A, after which all additions 
to A are obtained by applying a one-increment to the last nonzero or higher index 
in a multi-index already in A. 

An example of a relevant index set for which it is not straightforward to modify 
Algorithms and is the total degree with factorial correction TD-FC index set 
considered by Beck et al. [B]: 
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p = 2 

p = 4 

p = 6 

1/M 

M 

= 2 

.66 

.67 

.67 

.50 

M 

= 4 

.23 

.22 

.24 

.25 

M 

= 6 

.14 

.17 

.20 

.17 


p = 2 

p = 4 

p = 6 


M 

= 2 

9 

38 

no 


M 

= 4 

25 

255 

1519 


M 

= 6 

49 

924 

9324 



M 

= 2 M 

= 4 

M = 6 



K 800 60 20 





Number of chaos polynomials 


Figure 1. Timing results for a diffusion coefficient of the 
form (4.4) and an isoTD(M, iti) index set with different values of 
M,p, and K. Left: From top to bottom, estimated value of e in 
the complexity estimate 0 (|Wh||A|^+'^) together with the theoret¬ 
ical value the number of (sparse) matrix additions required 

for the computation of the summed matrices, i.e., I^(m)Ij 

and the maximum value of K considered when estimating e for dif¬ 
ferent M and p. Right: Time required to compute the neighbour 
matrices for M = 6 and p = 2,4,6 together with the line used to 
estimate e. The plot markers correspond to different values of K. 


Definition 4.1 (Total degree with factorial correction index set). Let N G N, 
pL G K^, and e > 0. The TD-FC index set is given by 

(4.5) 

TD-FC(7V, p, e) = |a e (N“)c | - log ^ < |'-logel,a„ = 0,n > ivj , 

where al = 

We cannot use the presented algorithms directly for the TD-FC index set as it 
is possible that if a ^ TD-FC(fV, p, e) then a + P G TD-FC(A^,/i, e), where /3 G 
(N)f’)c; for example, the multi-index (2,0,0,.. .) does not belong to the index set 
TD-FC(2, (1,1), but (2, 1, 0,0,.. .) does. 

To conclude this section, let us make a short comment about the general case. 
Assume we are given an index set, say A, with no topological information, i.e., 
the actual multi-indices in A have been generated, and we are asked to construct 
the neighbour matrices. Clearly Algorithm (or some variant of it) cannot be 
used directly as we do not have the tree structure provided by the parenthood 
relations available (or some other similar data structure). One option would be to 
construct a suitable tree structure, e.g., a KD-tree, and then use a modified version 
of Algorithm to construct the neighbour matrices. Constructing a KD-tree is a 
loglinear operation, and hence the complexity estimate for the construction of the 
neighbour matrices would not change in big O notation as 0(|A| log |A| -F |A|^+') = 
0(|A|^+'^). However, in practice it could be more efficient to first construct a 
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hash table for A using a suitable hash function and then use the hash table in 
Algorithmic to search elements in A. Employing hash functions in the construction 
of the stochastic moment matrices is in itself interesting and clearly out of the scope 
of this work, and hence we leave these type of considerations for future work. 


5. Numerical experiments 


We consider solving the following one-dimensional elliptic model problem in the 
Legendre case: find random field u S L^(r, iLg (D)) such that 


f -(a(y,x)u(y,a:)')'= 1 inD = (-l,l), 
1 u{y,-l) = u(y, 1) = 0 


holds for a given diffusion coefficient a(y, x) for which amin and Umax as in condition 
(1.2) exist. 

We consider the space-independent diffusion coefficient 


a(y,x) = a(y) = 1 -b 


M 


1 -I- ^ m ^yr, 


m—1 


and the space-dependent diffusion coefficient 


(5.2) 


a(y,a;) = 1 + 


M 


1 + m ® sin(m7ra;)i/„ 


with different values of M, s, and p as summarized in Table We have fixed 
the overall form of the diffusion coefficient and only varied the number of used 
random variables and the power in the diffusion coefficient; we could have consid¬ 
ered other types of diffusion coefficients, with similar conclusions. For the spatial 
discretization, we employ 20 p-type elements of order four. 

In the first example, we use a space-independent diffusion coefficient for which 
we have the exact solution available, and hence we can compare the performance of 
the various index sets against the exact result. The second and third examples are 
space-dependent for which we do not have a simple analytical solution available, and 
thus we rely on a case-specific “accurate enough” sGFEM solution as the “exact” 
solution. 

Even though the first example has a higher power in the diffusion coefficient 
than the two other examples, the small number of random variables and space- 
independency make it considerably easier to solve than the second or third example. 
Similarly, the second example is easier to solve than the third example as it has 
fewer random variables and a smaller power in the diffusion coefficient, which results 
in more sparse system matrices. In summary, we expect to require larger index sets. 



Diffusion coefficient type 

M 

s 

P 

Example 1 

space-independent 

2 

3/2 

6 

Example 2 

space-dependent 

4 

3/2 

2 

Example 3 

space-dependent 

6 

3/2 

4 


Table 1. Summary of the diffusion coefficient type and parame¬ 
ters M, s, and p in Examples 1-3. 
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i.e., more polynomials from the polynomial chaos, in the second and third examples 
than in the first one. 


5.1. Design of the numerical experiments. We denote Hq{D) with V and 
equip it with the gradient norm IMIu = l|Vu ||i 2 (£i),VR G V, which induces the same 
topology as the standard 7J-norm due to the homogeneous Dirichlet condition. The 


solution to (5.1) can be represented in terms of the Legendre polynomials as 


(5.3) 


i(y,a:)= ^ Uf,[x)L^{y) 

M6(N5°)c 


in the topology of Lp(r), where the Legendre “coefficient” G Hq{D) is given by 



u{y,x)Lf,{y)p{y)dy. 


By considering only a finite number of terms from ( |5.3[ ), we obtain an approximation 
for u 

u{y,x) « UA{y,x) = ^ u^{x)Lf,{y), 
msa 


where A is a finite multi-index set. The projection error is given by 


\\u-ua\\12^t,v) = XI 

which follows from Parseval’s equality and the completeness of in 

LliT). 

Based on na, it is assumed that for our numerical examples the following in¬ 
equality for the Legendre coefficient norms holds for s > 1 and any e > 0 when M 
approaches infinity: 


(5.4) 


El 

fi^A 




<C|A| 


1 —s+e 


CIAI 


1-S 


where (7 is a suitable constant. Inequality (5.4) implies the (weaker) inequality for 


(5.5) 




1 —S+e 


C|A| 


l-s 


We compare numerically observed Legendre coefficient convergence rates to (5.5) 
in Section [53] 


By solving (1.101 using a hnite multi-index set A (and a fine enough FEM grid), 
we obtain an sGFEM approximation for u 


u(y,x) « UAiy,x) = X +(a^)+(y), 

AtGA 

where {{ipIpgA are approximations for the Legendre coefficients We ob¬ 

tain the best M-terms approximation for u by first computing the (approximative) 
Legendre coefhcients of u in a sufficiently large index set A, i.e., {upjpgA, and then 
ordering the coefhcients in a non-increasing order based on their E-norm and tak¬ 
ing the first M-terms from the resulting ordered sequence. The corresponding M 
multi-indices form the bestM index set. 
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We report the performance of the TP, TD, and bestM index sets by considering 
the convergence of the relative errors of mean and variance (percentage), i.e. 


|E(ua) 


E(u) 


X 100 and 


Var({tA) — Var(r 


\L^D) 


\L^(D) 


Var('u)|| 


X 100, 


L^(D) 


when the number of multi-indices in the corresponding index set is increased. Here, 
“variance” is to be understood as the diagonal of the corresponding covariance 
function and the L^(Il)-norms are evaluated using built in numerical integration 
routines of Mathematica m- We denote the case-specific sGFEM solution func¬ 
tioning as the “reference” exact solution with u even though ua with a suitable 
index set A would be more appropriate. 

The weights for the dimensions used in the anisotropic index sets are computed 
numerically using ID analyses as in 011]. For each random variable 1 < m < M, 
we consider the subset = {/r S A : = 0 if i 7 ^ m, = 0 , 1 , 2 ,...}, where 

A C (Ng°)c is a sufficiently large index set. We assume the decay of the Legendre 
coefficients corresponding to A^ to be of the form ||M;i||v and we 

obtain the rate gm by linear interpolation on the quantities log ||m^||v. 

The numerical results for the experiments are obtained in the following way. 
First, the solutions corresponding to the isoTP(M, AT) and isoTD(M, AT) index 
sets, i.e., the isoTP(M, AT) and isoTD(M, AT) solutions, are computed by increasing 
AT from zero to some suitable maximum value. Then the weight vector g is ap¬ 
proximated from one of the computed sGFEM solutions as explained above. Once 
the weights are available, the anisotropic solutions are computed up to a suitably 
large index set. If the exact solution is not available, one of the computed sGEEM 
solutions is selected as the reference solution, and the relative errors in mean and 
variance are computed either against the exact or the reference solution. Einally, 
the best M-terms sGFEM solutions are computed as explained above by using the 
sGEEM solution corresponding to the largest used index set of the best performing 
index set type. 

In addition to the convergence plots of the relative errors in mean and variance, 
it turns out useful to visualize the norms of the Legendre coefficients of an sGEEM 
solution in descending order. From this kind of graph, it is possible to see how much 
excess multi-indices, i.e, multi-indices whose corresponding Legendre coefficients 
have small norms, are included in the index set. The “better” the index set type, 
the less excess multi-indices are included in the index set when the size of the index 
set is increased, i.e., the “tail” of the index set stays small. 

In the first example, where we consider only two random variables, we also illus¬ 
trate the largest considered index sets to give better understanding of the different 
index set types. Moreover, by plotting the norms of the Legendre coefficients as 
a function of the corresponding multi-indices, we are able to visualize why certain 
index sets are performing better than the rest of the index sets in that specific case. 


5.2. Space-independent case: Example 1. As the first example, we consider a 
space-independent diffusion coefficient with two random variables 


(5.6) 


iiy,x) = a{y) = 1 + 


-I 6 


1 + 


E 

m—l 


m 
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In this case, the exact solution to (5.1) is given by 

1-x^ 


and the main results are shown in Figures and 

Figure illustrates the convergence of the sGFEM solutions toward the exact 
mean (left) and variance (right) for the first experiment as the number of chaos 
polynomials, i.e., the number of multi-indices in the index set, is increased. From 
the figure, we see that the best results are obtained, as expected, with the bestM 
index set followed by the aTP, aTD, isoTD, and isoTP index sets. The anisotropic 
index sets are constructed using the weight vector g « (1.16,3.82), which is ap¬ 
proximated from the isoTD(2,19) solution, and the bestM index set is constructed 
based on the largest aTP solution. 

From Figure we see that the relative error in mean (variance) when using 
approximately 200 chaos polynomials is roughly 10“^% (10“^%) in the isoTP case, 
10“’^% (10“®%) in the aTP case, and 10“®% (10“^%) in the best M-terms case, and 
hence using different index set types may give dramatically better results. The 
aTP index set solutions follow first quite nicely the bestM index set solutions until 
the relative error between the index sets is about one decade at approximately 200 
multi-indices, indicating that the aTP index set is close to being optimal at least 
when the size of the index set is not too large. Recall, however, that the bestM 
index set is constructed based on the aTP solution with the largest index set. 

The solutions corresponding to the other index sets than to the isoTD index 
set converge relatively smoothly both in mean and variance whereas the solution 
corresponding to the isoTD index set is oscillating when considering the variance 
and stalling at every other increment when considering the mean. This reminds 
us that employing more computational resources to increase the size of the index 
set only slightly does not necessarily give better results as the solution may be 
oscillating even though it converges asymptotically. 

Figure shows the logarithm of the Legendre coefficient norms for the largest 
used isoTP index set, i.e., isoTP(2,24), on the left and the largest used index 
sets on the right for the first example. Notice that the range in the left image 
is different from the images on the right and that the left image and the top left 
image on the right correspond to each other; the latter shows the index set and 
the former portrays the log-norms of the corresponding Legendre coefficients. In 
a two-dimensional case, the isotropic/anisotropic tensor product and total degree 
index sets can be visualized as squares/rectangles and isosceles right triangles/right 
triangles, respectively. By looking at the log-norms shown in Figure we see that 
expanding the index set in the direction of the horizontal axis is more profitable 
than in the direction of the vertical axis. This agrees with the results in Figure 
the anisotropic index sets are performing better than the isotropic index sets. 

Figure 1^ illustrates the evolution of the ordered Legendre coefficient norms for 
the index sets considered in the first experiment as the dimension of the polynomial 
space increases. It seems that the best performing index sets have the smallest tails 
and that the tails grow slower for the anisotropic index sets than for the isotropic 
index sets, giving a new confirmation for the results in Figure]^ Notice that the 
largest index sets have almost equal relative errors as depicted by Figure and 
hence comparing the largest index sets is fair. If we want to decrease the size of 
an index set for which we have already computed the sGFEM solution, the ordered 




22 


HARRI HAKULA AND MATTI LEINONEN 


isoTP 
t— aTD 


isoTD ^ aTP 
f best M-terms 



Number of chaos polynomials 


isoTP 
■»— aTD 


isoTD ^ aTP 
^ best M-terms 



Number of chaos polynomials 


Figure 2. Convergence of the TP, TD, and bestM index set solu¬ 
tions toward the exact mean (left) and variance (right) for the first 
experiment as the dimension of the polynomial space increases. 


The diffusion coefficient is given by (5.6) and the anisotropic index 
sets are constructed using the weight vector g « (1.16, 3.82) which 
is approximated from the isoTD(2,19) solution. 
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Figure 3. The logarithm of the Legendre coefficient norms for 
the isoTP(2, 24) solution for the first example (left) and the largest 
used TP and TD index sets (right). The multi-indices in the index 
sets of the first example can be given as a = (ai, a 2 , 0,0,...), 
where Q!i,a 2 £ No- Here, the horizontal axis corresponds to ai 
and the vertical axis corresponds to a 2 - Notice that the range in 
the left image is different from those in the images on the right and 
that the left image and the top left image on the right correspond 
to each other. 


log-norm plot can be used to select a suitable threshold level for the multi-indices in 
the smaller index set. Finally, we have included an approximated convergence line 
with rate r = 2.27 in the log-log plots of Figure]^ to allow a comparison between 
the Legendre coefficients in the different examples; the larger the rate, the better 
approximations we are able to obtain with small index sets. The rate is obtained 
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isoTP 



10 ° 10 ' 10 ° 


isoTD 



10 ° 10 ' 10 ° 


aTP aTD 




10 ° 10 ' 10 ° 10 ° 10 ' 10 ° 


Figure 4. Evolution of the ordered Legendre coefficient norms 
(vertical axis) for the index sets considered in the first experi¬ 
ment as the dimension of the polynomial space increases (hori¬ 
zontal axis). We have also included an approximated convergence 
line with rate r = 2.27 in the images. 


from a least-squares fit to the 100 largest Legendre coefficients in the aTP case 
with the largest index set and it is computed in a similar way for the subsequent 
examples. Moreover, see Section 5.5 for a comparison of the observed Legendre 
coefficient convergence rates to the theoretical results. 


5.3. Space-dependent case: Example 2. In the second example, we use a space- 
dependent diffusion coefficient with four random variables 


(5.7) 


i(y,a:) = 1 + 


1 -I- m sm(rmTx)y„ 


m—l 


and we consider the sGFEM solution corresponding to the aTD index set with 4078 
chaos polynomials as the reference or “exact” solution as no simple form for the 
exact solution is available. 

Figuresandwhich are organized in the same way as Figures]^ andrespec¬ 
tively, for the first experiment, illustrate the convergence of the mean and variance 
of the sGFEM solutions toward the reference solution and the evolution of the 
ordered Legendre coefficients for the different types of index sets for the second 
example. Here, the weight vector g « (2.40,4.17,5.37,6.38) is approximated from 
the isoTD(4,17) solution, and the best M-terms index set is constructed based on 
the largest aTD solution. 

From Figure we see that the relative error in mean (variance) when using 
approximately 2 • 10^ chaos polynomials is roughly 10“®% (10“^%) in the isoTP 
case, 10“®% (10“^%) in the isoTD and aTP cases, 10“^®% (10“^°%) in the aTD 
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isoTP 
t— aTD 


isoTD ^ aTP 
f best M-terms 


isoTP 
■»— aTD 


isoTD ^ aTP 
^ best M-terais 




Number of chaos polynomials 


Number of chaos polynomials 


Figure 5. Convergence of the TP, TD, and bestM index set so¬ 
lutions toward the mean (left) and variance (right) correspond¬ 
ing to the aTD solution with 4078 chaos polynomials for the sec¬ 
ond experiment as the dimension of the corresponding polyno¬ 
mial space increases. The diffusion coefficient is given by (5.7) 
and the anisotropic index sets are constructed using the weight 
vector g « (2.40,4.17,5.37,6.38) which is approximated from the 
isoTD(4,17) solution. 


case, and 10“^^% (10“^^%) in the best M-terms case. As in the first example, the 
different index sets show dramatic differences in performance. Even though the 
aTP index set gave the best performance in the first example, the performance of 
the aTP index set is only on par with the generic isoTD index set in the second 
example. Finally, the isoTP index set performs poorly when compared with the 
other index sets. 

As in the first example, the best performing index sets have the slowest growing 
tails as illustrated by Figure This suggests the following heuristic method for 
selecting the index set to use in a practical application. First, use different index 
set types and compute the sGFEM solutions for these index sets using relatively 
small number of multi-indices. Finally, choose the index set which looks the most 
promising based on the evolution of the ordered Legendre coefficient norms. How¬ 
ever, this heuristic method is only feasible if the required sGFEM solutions can be 
computed in a reasonable time. 

The obtained convergence rate for the second example, r = 2.03, agrees with 
the convergence results of Figures and as the rate in the second example is 
lower than in the first example, we expect that larger index sets are required for a 
satisfactory converge of the solutions in the second example. See Section [57^ for the 
comparison of the observed Legendre coefficient convergence rates to the theoreticaf 
resuits. 
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isoTP 



10 ° 10 ' 10 ^ 10 ° 10 ^ 


isoTD 



10 ° 10 ' 10 ° 10 ° 


aTP 



10 ° 10 ' 10 ° 10 ° 


aTD 



10 ° 10 ' 10 ° 10 ° 


Figure 6. Evolution of the ordered Legendre coefficient norms 
(vertical axis) for the index sets considered in the second experi¬ 
ment as the dimension of the polynomial space increases (horizon¬ 
tal axis). We have also included an approximated convergence line 
with rate r = 2.03 in the images. 


5.4. Space-dependent case: Example 3. As the third and final example, we 
use a space-dependent diffusion coefficient with six random variables 


(5.8) 


a(y,ai) = 1 -P 


1 -P m sin(m7ra;)y„ 


and we consider the sGFEM solution corresponding to the aTD index set with 5909 
chaos polynomials as the reference or “exact” solution as no simple form for the 
exact solution is available. The weight vector g « (1.52,3.56,4.95,5.76,6.4,6.85), 
used in the anisotropic solutions, is approximated from the isoTD(6,9) solution. 

Figure depicts the convergence of the relative errors for mean and variance 
in a similar way as in the previous two examples. It seems that the aTD index 
set clearly outperforms the aTP, isoTD, and isoTP index sets, the aTP index set 
exhibits slightly better performance than the isoTD index set, and the isoTP index 
set, as in the previous two examples, gives the worst performance. In this example, 
we were limited by our computing resources to around 5 • 10^ multi-indices due 
to more full system matrices caused by the higher power and number of random 
variables in the diffusion coefficient than in the second example. Hence, the largest 
index sets do not correspond to similar errors but the cardinalities of the largest 
index sets are about the same. 

Figure which shows the ordered Legendre coefficient norms for the third ex¬ 
ample, agrees with the previous two examples: the slower the tail is growing for an 
index set type, the better it is performing. As we have increased both the number 
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isoTP 
1— aTD 


isoTD — ^— aTP 

r best M-terms 



10 ' 10 ^ 10 ’ 
Number of chaos polynomials 


isoTP 

aTD 


isoTD ^ aTP 
^ best M-terms 



10 " 10 ' 10 " 10 ’ 
Number of chaos polynomials 


Figure 7. Convergence of the TP, TD, and bestM index set solu¬ 
tions toward the mean (left) and variance (right) corresponding 
to the aTD solution with 5909 chaos polynomials for the third 
experiment as the dimension of the corresponding polynomial 
space increases. The diffusion coefficient is given by (5.8) and 


the anisotropic index sets are constructed using the weight vector 
g « (1.52, 3.56,4.95, 5.76, 6.4,6.85) which is approximated from the 
isoTD(6,9) solution. 


of random variables and the power in the diffusion coefficient, it is natural that the 
approximated convergence line in the log-log plots of Figure has a smaller rate, 
r = 1.43, than in the previous two examples: more chaos polynomials are required 
to take into account the increased number of random variables and cross terms in 
the diffusion coefficient. According to Figurethe isoTP index set and the isoTD 
and aTP index sets start to perform worse than the aTD index set at around 10^ 
and 5 • 10^ multi-indices, respectively, which is in line with the results in Figure 


5.5. Legendre coefficient convergence rates. To conclude the numerical exper¬ 
iments, we compare the observed Legendre coefficient convergence rates to theoret¬ 
ical results, i.e., we examine how close to the theoretical rate = s — 1 from (5.5) 
the Legendre coefficients converge when M is large enough; recall that in our previ¬ 
ous three experiments s = 3/2. We have listed approximated convergence rates for 
space-dependent diffusion coefficients of the form (5.2) with different values of M, s, 
and p in Table|^ Considerably larger values of M than shown on Tablej^could have 
been considered for the case p = 1. However, as this would not have been feasible 
for p = 2, 3, we restricted M below 9. The rate for particular values of M, s, and 
p in Table was obtained by first computing the solution using the isoTD(M, r) 
index set, where r is the smallest value for which | isoTD(M, r)| > 3 • 10^, after 
which the rate is computed as in the previous examples using Legendre coefficients 
from the 10th largest up to the 100th largest. The largest Legendre coefficients 
were excluded for stability reasons. 

The values in Table agree relatively well with the related theory: the rates 
corresponding to s = 6, 8 with M = 8 and p = 1,2,3 are close to the expected 
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10 " 10 ' 10 ^ 10 ^ 10 “ 10 ' 10 ^ 10 “ 




10 " 10 ' 10 “ 10 “ 10 " 10 ' 10 “ 10 “ 

Figure 8. Evolution of the ordered Legendre coefficient norms 
(vertical axis) for the index sets considered in the third experi¬ 
ment as the dimension of the polynomial space increases (horizon¬ 
tal axis). We have also included an approximated convergence line 
with rate r = 1.43 in the images. 


values of 5 and 7, respectively. For the cases s = 2,4, the obtained rates of ap¬ 
proximately r 2 = 2 and r 4 = 3.5 are slightly better than the expected values of 1 
and 3 predicted by |14j . However, by increasing M the rates keep approaching the 
predicted values: for example, M = 50 gives rates r 2 = 1-66 and r 4 = 3.26 (not 
shown), respectively, and hence we expect that by further increasing M the rates 
would eventually approach 1 and 3, respectively. In summary, the larger s is, the 
faster we expect the converge rate to approach s — 1 when M is increased. 

According to Table the power p does not essentially affect the rate. This is ex¬ 
pected: if we expanded the diffusion coefficient, the terms resulting from increasing 
the value of p have smaller overall effect on the diffusion coefficient. However, as 
the estimated rate depends on the actual Legendre coefficients used in the compu¬ 
tation of the rate, we have plotted the ordered Legendre coefficients corresponding 
to the last row of Table i.e., M = 8 to confirm the results: Figure shows the 
corresponding four sets (s = 2,4, 6, 8) of almost overlapping curves {p = 1, 2, 3). As 
the ordered Legendre coefficients have the same profile for each fixed value of s, we 
conclude that the power p does not (essentially) affect the asymptotic rate. 

By ignoring the power p and the type of the diffusion coefficient, we know from 
Sections |5.2f|5.4| that the values M = 2,4,6 with s = 3/2 correspond to the rates 
^ 3/2 = 2.27, 2.03, and 1.43, respectively. These numbers are in line with the results 
presented in Tablej^for the case s = 2. By increasing M the rate keeps decreasing: 
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p = 3 

M s=2 4 6 8 


p = I p = 2 


M 

s=2 

4 

6 

8 

M 

s=2 

4 

6 

8 

2 

7.2 

9.3 

11.0 

12.5 

2 

6.1 

8.0 

9.6 

11.0 

4 

2.9 

4.5 

5.8 

7.2 

4 

2.6 

4.0 

5.4 

6.8 

6 

2.2 

3.7 

5.2 

6.9 

6 

2.1 

3.5 

5.1 

6.8 

8 

2.0 

3.5 

5.1 

6.8 

8 

1.9 

3.5 

5.1 

6.8 


5.1 

6.8 

8.4 

2.3 

3.7 

5.1 

1.9 

3.4 

5.1 

1.9 

3.4 

5.0 


Table 2. Legendre coefficient convergence rates for the space- 
dependent diffusion coefficient (5.2) with different values of M, s, 
and p. With M = 8 and p = 1,2,3 the expanded form of the 
diffusion coefficient has 9, 45, and 165 terms, respectively. 


9.6 

6.6 
6.8 
6.7 



10 " 10 ' 10 " 10 " 


Figure 9. The ordered Legendre coefficient norms (vertical axis) 
corresponding to the last row of Table i.e., M = 8. The four 
sets of three almost overlapping lines correspond to the four and 
three different values of s and p, respectively, in Table 


for example, M = 50 gives rate r ^/2 = 1.29 (not shown), and hence we expect that 
by further increasing M the rate would eventually approach 1/2. 


6. Discussion 


In this paper, we have presented algorithms for constructing stochastic moment 
matrices of the form (ol) in an efficient manner for certain index set types allowing 
the modeler to consider non-affine diffusion coefficients of the form (1.4). 

Algorithms and can be used (with possible slight modifications) to the effi¬ 
cient construction of hierarchical index sets together with the neighbour matrices 
required in the stochastic moment matrix generation. The algorithms can be used 
with index sets which can be constructed by adding an initial multi-index to the 
index set, after which all additions to the index set are obtained by applying a cer¬ 
tain kind of one-increment to a multi-index already in the index set. Most common 
index sets can be constructed in this way, with the exception of, for example, the 
factorial corrected total degree (TD-FC) index set. We leave the generalization of 
Algorithms [^and|^to index set types where any one-increment is allowed for future 
studies. 
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In the numerical examples, we reported the performance of various index set 
types for one space-independent and two space-dependent diffusion coefficients in 
a one-dimensional spatial setting. Moreover, we compared the observed Legendre 
coefficient convergence rates to theoretical results: the observed Legendre coefficient 
rates agreed well with the theoretical results in [14]. 

Using the isotropic total degree index set to compute the dimension weights for 
the anisotropic index sets as in HIT] worked well; the anisotropic index sets outper¬ 
formed the isotropic index sets in all considered examples. In the first experiment, 
the anisotropic tensor product index set gave the best results, whereas in the sec¬ 
ond and third experiment the anisotropic total degree index set performed the best. 
Hence, it might be difficult to predict a priori which index set type would be the best 
choice without a possibly extensive analysis of the problem in question. According 
to our numerical experiments, a promising heuristic method to select the index 
set type is to first compute sGFEM solutions for various index set types and then 
select the one with the slowest growing “tail” of Legendre coefficient norms. Other 
possibility is to compare the sGFEM solutions, in some suitable sense, with a set of, 
for example, FEM solutions corresponding to different realizations of the diffusion 
coefficient and then select the best performing index set type. The (numerical) 
consideration of these ideas is left for future studies. 


Appendix A. Legendre and Hermite polynomials 


When constructing the stochastic system matrices for (1.6) or (1.10) in the set¬ 
ting of generalized polynomial chaos (gPG) expansions, see, e.g., we often 
encounter products of orthogonal hypergeometric polynomials. It turns out useful 
to linearize these products using a suitable polynomial family, i.e., solve the (gen¬ 
eralized) linearization problem for the polynomial families involved. Moreover, it 
is also useful to solve the inversion problem for each of the polynomial families, 
i.e., express where A: is a nonnegative integer, using the polynomial family in 
question; see [3| and the references therein for more information about linearizing 
expressions containing (hypergeometric) polynomials. 

In this paper, we consider only the Legendre and Hermite polynomial families as 
those are typical choices in the sGFEM setting. Hence, we consider the lineariza¬ 
tion and inversion problem only for these (hypergeometric) orthogonal polynomial 
families. However, similar results to the ones shown below can be obtained for 
other (hypergeometric) orthogonal polynomial families [3|. 

In the next two sections, we recall the linearization, triple integral, and inversion 
formulas for the Legendre and Hermite polynomials. We note that even though the 
presented results are classical, the actual formulas are rarely seen in the sGFEM 
literature. 


A. I. Legendre polynomials. The orthonormal univariate Legendre polynomials 
are given explicitly by the Rodrigues’ formula: 


Definition A.l (Legendre polynomials). Let m G Nq. The mth orthonormal 
univariate Legendre polynomial is given by 


■\/2m”+T d” 




(y) 


2^ ml 
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with the orthonormality relation 


J Lk{y)Li{y)-dy = Ski, 
where Ski is the Kronecker’s delta. 

Theorem A.2 (Linearization formula for the Legendre polynomials). The product 
of two univariate Legendre polynomials Lk{y) and Li{y) can be expanded as a sum 
of Legendre polynomials: 

k+l 

^k{y)Ll(^y) = ^ ^ -^{23 is even}' 

(A.l) m=|fc-Z| 


A{g - k)A{g - l)A{g - m) 


A{g) 

(2k + 1)(2Z + 1) (2to + 1) 


25 + 1 


Lmiy), 


where g = {k + I + m)/2 and A{n) = /2". 

Proof. See [T]. 


□ 


The coefficient in the series (A.l I can be written in a simpler form by employ¬ 
ing the 3jm-symbol. We need only the specific case with zero magnetic quantum 
numbers; see m for more information about the 3jm-symbols. 

Definition A.3 (3jm-Symbol). The 3jm-symbol with zero magnetic quantum num¬ 
bers and k,l,m G Nq is defined as 

2 /n . F _l n2 


k I m 
0 0 0 


(25 - 2fc)!(25 - 20!(25 - 2 to)! 

' (25 + 1)! ~ 


9'- 


{g - k)\{g - iy.{g - m)'.\ 


if 2g is even and \k — l\ < m < k + l; otherwise (g o o*zero. Here, g is as in 
Theorem IA.2I 


By noticing that 


^{2g is even}l-{|fc —/|<m<fc-t-/} 


A{g-k)A{g-l)A[g-m) 1 


^(5) 


25+1 


k I m 
0 0 0 


we can write (A.l) as 

(A.2) 

where 


k-\-l 


Lk{y)Li{y)= X! L{k,l,m)Lm{y), 

m—\k—l\ 


L{k,l,m) = \J {2k + 1){21 + l)(2m + 1) 


k I m 
0 0 0 


Using the linearization formula ( |A.2 ), we can compute the triple product integral 
of the Legendre polynomials. 

Theorem A.4 (Triple product integral of the Legendre polynomials). The triple 
product integral of univariate Legendre polynomials Lk{y), L^y), and Lm{y) is given 
by 

/ Lk{y)Li{y)Lm{y)p{y)dy = L{k, I, to), 
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where 


Proof. 


L{k,l,m) = \/ {2k + 1){21 + l)(2m + 1) 


k I m 
0 0 0 


/ I k+l 

Lk{y)Li{y)Lm{y)p{y)dy = ^ L{k, I, n)Snm = L{k, I, m), 

^ n—\k—l\ 

where we have first used ( A.2[ ) and the orthonormality of the Legendre polynomials 
and then the selection rule \k — l\ < m < k + I ot L{k, l,m). □ 

Finally, we give the inversion formula for the univariate Legendre polynomials. 

Theorem A.5 (Inversion formula for the Legendre polynomials). The solution to 
the inversion problem for the univariate Legendre polynomials is given by 

k 


y 


= J2lnLn{y), 


n—0 


where 


IZ = 1 


f^\ I m --(fc —n—1)!! „ , 

{k-n IS even} ^ 1) 11 ’ 0 < U < fc, 

and !! denotes the double factorial. 

Proof. The result can be obtained, e.g., by using the techniques in [3]. 


□ 


A.2. Hermite polynomials. We start by defining the Hermite polynomials and 
then give the linearization, triple product integral, and inversion formulas for them. 

Definition A.6 (Hermite polynomials). Let m G Nq. The ruth orthonormal (prob- 
abilists’) univariate Hermite polynomial is given by 


Hm{y) = --^exp{y^/2)—exp{-y^/2), 
Vm! dy"* 

with the orthonormality relation 


TO = 0,1,... , 


Jr v27r 

Theorem A.7 (Linearization formula for the Hermite polynomials). The series 
for the product of two univariate Hermite polynomials Hk{y) and Hi{y) is given by 

k-\-l 

Hk{y)Hi{y) = ^ H{k,l,m)Hm{y), 

m—\k — l\ 

where 

1 

r / ^ \ / I 

Hl^kjl^nT) is even} ^{\k — l\<m<k+l} 

with g = {k + I + to)/2. 

Proof. See [3]. □ 

By using Theorem |A.7[ we obtain the triple product integral of the Hermite 
polynomials in a similar way as in the Legendre case. 


g — mJ \g — mJ \g — k J 
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Theorem A.8 (Triple product integral of the Hermite polynomials). The triple 


product integral of univariate Hermite polynomials Hk(y), Hi{y), and is 

given hy 



JK. 

We conclude with the solution to the inversion problem for the Hermite polyno¬ 
mials. 

Theorem A.9 (Inversion formula for the Hermite polynomials). 

k 


y 






where 



0 < n < k. 


Proof. See [3]. 


□ 
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Appendix B. Pseudocode for generating the ts index set 


Algorithm 3: Pseudocode for generating the index set TS(^,c) and the corre¬ 
sponding parenthood relations PtC. 

Data: Vector p. S cq and real number e G (0,1). 

Result: Multi-indices TS(^,e) and parenthood relations PtC. 

1 initialize a to empty vector of integer pairs; 

2 push pair (a, 1) to the end of TS(/r, e); 

3 initialize Sa, m, p, and c to one; 

4 initialize n to two; 

5 create empty stacks for a, to, Cq,, and p; 


6 while True do 


7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 


if EaPm > £ then 

push a, TO, and £„ to the corresponding stack; 
push c to p-stack; 
increase to by one; 
continue; 

end 

if m-stack is empty then return TS(/i,£) and PtC; 

set a, TO, Ea, and p to the top value in the corresponding stack; 

pop a, TO, Ea, and p-stacks; 

if # of pairs in a = 0 OR a[—1][1] ^ to then push integer pair (to, 0) to 

the end of a ; 

increase a[—1][2] by one; 

multiply Ea with 

set c to n; 

increase n by one; 

push pair {a,Ea) to the end of TS(p,£); 

while # of elements in PtC < p do push an empty integer vector to the 
end of PtC; 

push c to the front of PtC[p]; 


24 end 


Appendix C. Example run of Algorithm [3] 

An example run of Algorithm with intermediate steps in the spirit of Algo¬ 
rithm [l] is shown in Figure [T0| and in Table for 

p= (2-\3-\4-\50-\60-\70-\...) and £ = 20"^ 

Figure [To| shows both the resulting index set TS(p, £) and the corresponding PtC 
information, and Table [^contains the intermediate steps, i.e., for each pass through 
the main loop in the algorithm, the most recent addition to the index set, the stack 
just before the next addition to the index set, and the considered child which was 
not added to the stack. 
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Figure 10. Index set TS(^,e) and the corresponding PtC information for /i = (2“^, 3“^, 4“^, 50“^, 60“^, 70“^,...) 
and e = 20“^. The vertices correspond to the multi-indices in TS(/i,e). The vertex labels show the generation order 
(top), the first three indices (the rest are zeros) (middle), and the (rounded) value (bottom) of the multi-indices. 
The directed edges correspond to the parenthood relations between the multi-indices, i.e., they display the information 
stored in PtC. The edge labels indicate the order of the children of, say, a in the vector PtC[a]. 


KARRI HAKULA AND MATTI LEINONEN 
































# 

Addition 

Stack 

Considered 

# 

Addition 

Stack 

Considered 

1 

( 0 , 0 , 0 ) 

1 

( 1 , 0 , 0 ),( 0 , 1 , 0 ), ( 0 , 0 , 1 ) 

( 0 , 0 , 0 , 1 ) 

0.02 

9 

( 1 , 1 , 0 ) 

0.17 

( 2 , 0 , 0 ),( 1 , 2 , 0 ) 

( 1 . 1 , 1 . 0 ) 

0.04 

2 

( 0 , 0 , 1 ) 

0.25 

( 1 , 0 , 0 ),( 0 , 1 , 0 ), ( 0 , 0 , 2 ) 

( 0 , 0 , 1 , 1 ) 

0.01 

10 

( 1 , 2 , 0 ) 

0.06 

( 2 , 0 , 0 ) 

( 1 , 3 , 0 , 0 ) 

0.02 

3 

( 0 , 0 , 2 ) 

0.06 

( 1 , 0 , 0 ),( 0 , 1 , 0 ) 

( 0 , 0 , 3 , 0 ) 

0.02 

11 

( 2 , 0 , 0 ) 

0.25 

( 3 , 0 , 0 ),( 2 , 1 , 0 ), ( 2 , 0 , 1 ) 

( 2 , 0 , 0 , 1 ) 

0.01 

4 

( 0 , 1 , 0 ) 

0.33 

( 1 , 0 , 0 ),( 0 , 2 , 0 ),( 0 , 1 , 1 ) 

( 0 , 1 , 0 , 1 ) 

0.01 

12 

( 2 , 0 , 1 ) 

0.06 

( 3 , 0 , 0 ),( 2 , 1 , 0 ) 

( 2 , 0 , 2 , 0 ) 

0.02 

5 

( 0 , 1 , 1 ) 

0.08 

( 1 , 0 , 0 ),( 0 , 2 , 0 ) 

( 0 , 1 , 2 , 0 ) 

0.02 

13 

( 2 , 1 , 0 ) 

0.08 

( 3 , 0 , 0 ) 

( 2 , 2 , 0 , 0 ) 

0.03 

6 

( 0 , 2 , 0 ) 

0.11 

( 1 , 0 , 0 ) 

( 0 , 3 , 0 , 0 ) 

0.04 

14 

( 3 , 0 , 0 ) 

0.13 

( 4 , 0 , 0 ) 

( 3 , 1 , 0 , 0 ) 

0.04 

7 

( 1 , 0 , 0 ) 

0.5 

( 2 , 0 , 0 ),( 1 , 1 , 0 ),( 1 , 0 , 1 ) 

( 1 , 0 , 0 , 1 ) 

0.01 

15 

( 4 , 0 , 0 ) 

0.06 


( 5 , 0 , 0 , 0 ) 

0.03 

8 

( 1 , 0 , 1 ) 

0.13 

( 2 , 0 , 0 ),( 1 , 1 , 0 ) 

( 1 , 0 , 2 , 0 ) 

0.03 



Table 3. Intermediate steps in Algorithm when constructing the index set TS(/r, e) for = 

(2“^, 3“^, 4“^, 50“^, 60“^, 70“^,...) and e = 20“^. From left to right, the columns consists of the number of times we 
have passed through the start of the main loop, the first three indices (the rest are zeros) of the most recently added 
multi-index to TS(^,£:) with the corresponding (rounded) value /i“ shown below, multi-indices in the stack just before 
the next addition to the index set, and the considered child which was not added to the stack. 


W 

CJi 
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Appendix D. Pseudocode for generating the neighbour matrices 

Algorithm 4: Pseudocode for generating the neighbour matrices 
corresponding to the index set TS(/x, e). 

Data: Weights WA;, index set TS(^,£), and the corresponding vector real 
number £, and parenthood relations PtC. 

Result: Neighbour matrices {N'^}uj^Ws- 

1 for m = 1 to # of elements in MA do 

2 set w to Wh[to]; 

3 set Sw to /X™; 

4 for n = 1 to # of elements in TS(/x,e) do 

5 set ?7 to TS(^, e)[n][l]; 

6 set to TS(/x, £)[n][2]; 

7 set 7 to ?7 — re; 

8 if 7 [r][2] < 0 for some r OR e^few < e then continue; 

9 initialize k and u to one; 

10 for each pair p in ^ in the increasing position order do 

11 if p[2] = 0 then continue; 

12 initialize j to p[2]; 

13 if p[l] > u then 

14 decrease j by one; 

15 set k to PtC[fc][p[l] — w + 1]; 

16 setMtop[l]; 

17 end 

18 while j > 0 do 

19 set k to PtC[fc][l]; 

20 decrease j by one; 

21 end 

22 end 

23 set [a™] ^ ^ and [A™] „ to one; 

24 end 

25 end 

26 return 
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Appendix E. Modifications for the isoTD index set 

Modifications to Algorithmfor the isoTD index set. 

Data: Positive integer N and nonnegative integer K. 

2 push pair (a, 0) to the end of TD(A, K); 

3 initialize Cq, to zero; initialize m, p, and c to one; 

7 if Ea < K AND m < N then ... ; 

18 increase £„ by one; 

Modifications to Algorithmfor the isoTD index set. 

Data: Weights VIA, index set isoTD(Af, AT), and the corresponding K, and 
PtC. 

3 set Syj to y 1 

8 if length(i(;) > N OR 7[r][2] < 0 for some r OR Sr^ — > K then continue; 


Appendix F. Modifications for the aTP index set 


Modifications to Algorithmfor the aTP index set. 


Data: Positive integer N, nonnegative integer A, and weight vector p S 

2 push pair (a,0) to the end of aTP(A, A,/i); 

3 initialize Sa to zero; initialize m, p, and c to one; 

7 if Ea < K AND m < N then ...; 

10 increase m by one; set Ea to zero; 

13 while True do 


13 

18 


if m-stack is empty then return aTP(A, A, and PtC; 


increase by 


if £„ < A then break; 


18 end 


Notice that the lines 13 and 18 in Algorithm are replaced with multiple lines. 
Modifications to Algorithm for the aTP index set 
Data: Weights WA, multi-indices aTP(A, A,/i), and the corresponding N, A, 
/i, and PtC. 

3 remove this line; 

8 if length(w) > N OR 7[r][2] < 0 OR ^b[blPl'y[r][2] > A for some r then 
continue; 

15 set k to PtC[fc][—1 — aTP(A, A,^)[PtC[A:][—1]][1][—1][1] +p[l]]; 


Here, PtC[a][— m] denotes the mth element of the vector PtC[a] counting from 
the end. In the aTP case it might happen that even though the multi-indices 
(oi,..., Om, 0, 0,...) and (oi, Ij 0,...) belong to the aTP index set the 

multi-index (ai,...,am + 1,0,0,...) does not. Due to this fact, we had to in¬ 
troduce a new while loop between lines 13 and 18 of Algorithm to ensure that 
only valid multi-indices are added to the index set. Moreover, we had to process 
the position-value pairs in a different manner than in the other cases resulting in 
the difficult looking line 15 in Algorithm]^ On that line, the one step down the 
tree to the child with the length m is done from the end and not from the front. 
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If the dimension weights /x are given using the vector g, we need to replace 

on line 18 of the aTP modified Algorithm with gm and on line 8 of the aTP 

modified Algorithmwith and we consider AT as a nonnegative real number. 
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